---
title: "Replication of Vogel and Willems"
author: "Paw Hansen"
date: "`r Sys.Date()`"
output: html_document
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = F, message = F, warning = F, fig.height = 5, fig.width = 5, fig.align = "center")

packages <- c("here","readr","tidyverse", "broom", "cowplot",
              "psych", "forestmangr", "scales", "ggrepel", "janitor", "MASS")
library(tidyverse)

map(packages, library, character.only = TRUE)
theme_set(theme_cowplot())
```

# Import and recode data
```{r}
data <- read_csv(here("Data", "data_succes.csv")) %>% 
                   as_tibble() 

# Codes for recoding new variables
likert_turnover<- c("Helt uenig"=1,
                    "Unig"=2, 
                    "Hverken enig eller uenig"=3,
                    "Enig"=4, 
                    "Helt enig"=5)

likert_turnover_reverse <- c("Helt uenig"=5, 
                          "Unig"=4, 
                          "Hverken enig eller uenig"=3, 
                          "Enig"=2, 
                          "Helt enig"=1)

age_codes <- c("18-24 år"=1,
               "25-29 år"=1,
               "30-34 år"=2,
               "35-39 år"=2,
               "40-44 år"=3,
               "45-49 år"=3,
               "50-54 år"=4,
               "55-59 år"=4,
               "60-64 år"=5,
               "Mere end 64 år"=5)

male_codes <- c("Kvinde" = 0, "Mand" = 1, "Andet/ønsker ikke at svare" = NA)

treat_codes <- c("success_1"= "Recent_user",
                 "success_2" = "Recent_soc",
                 "success_3" = "Special_user",
                 "success_4"="Special_soc",
                 "success_5"="Control"
                 )

jobsat_codes <- c("1: Meget utilfreds"=1, 
                  "2"=2,
                  "3"=3,
                  "4"=4,
                  "5"=5,
                  "6"=6,
                  "7"=7,
                  "8"=8,
                  "9"=9,
                  "10: Meget tilfreds"=10)

data <- data %>% 
  mutate(
    Age_group=recode(age, !!!age_codes),
    Treatment =recode(FL_12_DO, !!!treat_codes),
    treatment_diko = case_when(
      Treatment == "Recent_soc" | Treatment == "Special_soc" ~ "Societal",
      Treatment == "Recent_user" | Treatment == "Special_user" ~ "Prosocial",
      Treatment == "Control" ~ "Control"),
    Male = recode(sex, !!!male_codes), 
    Jobsatisfaction =recode(jobsatisfaction, !!!jobsat_codes),
    Turnover_1 =recode(`turnover-intention_1`,!!!likert_turnover_reverse),
    Turnover_2 =recode(`turnover-intention_2`,!!!likert_turnover),
    Turnover_3 =recode(`turnover-intention_3`,!!!likert_turnover),
    duration = `Duration (in seconds)`
    ) %>% 
  dplyr::select(
         -contains("intention"), 
         -jobsatisfaction,
         -age,
         -sex,
         -FL_12_DO) %>% 
  clean_names() 
```

Turnover scale
```{r, results='hide'}
data <- data %>%
  mutate(turnover = dplyr::select(., turnover_1:turnover_3) %>%
           rowSums(na.rm = TRUE)) %>% 
  relocate(c("treatment", "jobsatisfaction", "turnover"))

data %>%
  dplyr::select(turnover_1:turnover_3) %>%
  psych::alpha(check.keys=TRUE) #Std. Cronbachs alpha user orientation scale: 0.89
```

```{r}
# Rescale
data <- data %>% 
  mutate(across(c(jobsatisfaction, turnover), scales::rescale, to = c(1, 7))
         )

# NA on key variables? 
data %>% 
  summarise(across(c(treatment, treatment_diko, jobsatisfaction, turnover_1), is.na)) %>% 
  map_dbl(mean)
```

# Findings

## Descriptives
```{r}
library(vtable)

descriptives <- data %>% 
  drop_na(treatment_diko, jobsatisfaction, turnover) %>% 
  dplyr::select(treatment_diko, male, age_group, education, edd_socialworker, workexperience_casew, tenure_curren_job) %>% 
  mutate(education = as.factor(education),
         age = as.factor(age_group))

descriptives_table <- sumtable(descriptives, 
         digits = 2, 
         fixed.digits = TRUE,
         out='return',
         group = 'treatment_diko') 

write.csv(descriptives_table, 
            file = here("Figures", "descriptives_table.csv"),
            quote = FALSE, 
            row.names = F)

# For counts (n)
data %>% 
  drop_na(treatment_diko, jobsatisfaction, turnover) %>% 
  count(treatment_diko) 

```

## Average Treatment Effects
```{r}
fit_jobsat <- lm(jobsatisfaction ~ treatment_diko, data = data)
```

```{r}
fit_turnover <- lm(turnover ~ treatment_diko, data = data) 
```

## Plots 

Jobsatisfaction for each treatment group, with means and +/- one standard deviation point ranges.
```{r, fig.show='hide'}
(summary_jobsat <- data %>% 
  drop_na(treatment_diko) %>% 
  group_by(treatment_diko) %>% 
  summarise(model = list(tidy(lm(jobsatisfaction ~ 1), conf.int = T))) %>% 
  unnest(model) %>% 
  rename(jobsatisfaction = estimate))

set.seed(1234)

(p_jobsat <- data %>% 
  drop_na(treatment_diko, jobsatisfaction) %>% 
  ggplot(aes(x = treatment_diko, y = jobsatisfaction)) + 
  geom_jitter(alpha = 0.5, width = 0.2, color = "steelblue") +
  geom_hline(yintercept = mean(data$jobsatisfaction[data$treatment_diko=="Control"], na.rm=T), linetype = "dotted") +
  geom_pointrange(data = summary_jobsat, aes(ymin = conf.low, 
                                             ymax = conf.high)) + 
  scale_y_continuous(breaks = seq(1, 7, 1)) + 
  scale_x_discrete(labels = c("Control\n (n = 84)", "Prosocial\n (n = 163)", "Societal\n (n = 165)")) + 
  labs(x = NULL, 
       y = "Jobsatisfaction")  
  )
```

```{r, fig.show='hide'}
(summary_turnover <- data %>% 
  drop_na(treatment_diko) %>% 
  group_by(treatment_diko) %>% 
  summarise(model = list(tidy(lm(turnover ~ 1), conf.int = T))) %>% 
  unnest(model) %>% 
  rename(turnover = estimate))

set.seed(1234)

(p_turnover <- data %>% 
  drop_na(treatment_diko, turnover) %>% 
  ggplot(aes(x = treatment_diko, y = turnover)) + 
  geom_jitter(alpha = 0.5, width = 0.2, color = "steelblue") +
  geom_hline(yintercept = mean(data$turnover[data$treatment_diko=="Control"], na.rm=T), linetype = "dotted") +
  geom_pointrange(data = summary_turnover, aes(ymin = conf.low, 
                                             ymax = conf.high)) + 
  scale_y_continuous(breaks = seq(1, 7, 1)) + 
  scale_x_discrete(labels = c("Control\n (n = 84)", "Prosocial\n (n = 163)", "Societal\n (n = 165)")) + 
  labs(x = NULL, 
       y = "Turnover intention"))

```

```{r}
library(patchwork)
(fig_1 <- p_jobsat / p_turnover + plot_annotation(tag_levels = 'A'))

ggsave(filename= "treatment_fig.tiff", 
       path="Figures",
       width = 8, 
       height = 8, 
       device='tiff', 
       dpi=700)
```


## As tables 
```{r}
nice_table <- function(model) {
  tidy(model, conf.int = T) %>% 
  round_df(digits = 2) %>% 
  mutate(term = str_remove(term, "treatment_diko"),
         confidence_interval = paste0("[", conf.low, " - ", conf.high, "]")) %>% 
  dplyr::select(-statistic, -conf.low, -conf.high)
}

model_jobsat <- nice_table(fit_jobsat)

write.csv(model_jobsat, 
            file = here("Figures","model_jobsat.txt"), 
            quote = FALSE, 
            row.names = F)

model_turnover <- nice_table(fit_turnover)

write.table(model_turnover, 
            file = here("Figures","model_turnover.txt"), 
            sep = ",", 
            quote = FALSE, 
            row.names = F)
```

# Manipulation check

Time spent completing the survey (t-tests)?
```{r}
data %>% 
  drop_na(treatment_diko) %>% 
  group_by(treatment_diko) %>% 
  filter(duration < 3600) %>% # Keep only if spent less than 60 minutes
  summarise(duration = list(tidy(t.test(duration/60)))) %>% 
  unnest(duration) %>% 
  round_df(digits = 2)
```

Filled out open text field with reflection?
```{r}
data %>% 
  drop_na(treatment_diko) %>% 
  group_by(treatment_diko) %>% 
  summarise(wrote_something = 1-mean(is.na(succes_essay))) %>% 
  round_df(digits = 2)

data %>% 
  drop_na(treatment_diko) %>% 
  group_by(treatment_diko) %>% 
  summarise(avg_n_characters = mean(str_length(succes_essay), na.rm=T))
```

# Probabilistic analysis 
```{r}
# Simulate estimation uncertainty
fit_turnover %>% 
  tidy()

set.seed(1234)
sim_eff_turnover <- mvrnorm(10000, mu = coef(fit_turnover), Sigma = vcov(fit_turnover)) %>%
  as_tibble()

# User defined function 
less_than <- function(variable, threshold) {
  mean(variable <= threshold)
}
```

## Prosocial impact
```{r}
# Sequence of the treatment effects sizes of interest
seq_treat_size <- seq(to = 1.1, from = -1.1, by = .001)

# Probability of getting each of those treatment effect sizes
probs <- seq_treat_size %>% 
  map_dbl(.f = less_than,
          variable = sim_eff_turnover$treatment_dikoProsocial)

# Bind them together for plotting
results <- cbind(seq_treat_size, probs) %>% 
  as_tibble()

vogel_willems_results <- 
  tribble(~study, ~ate,
          "Study 1", -0.63,
          "Study 2", -0.16,
          "Study 3", -0.42)

(p_prob_prosoc <- results %>% 
  ggplot(aes(x = seq_treat_size, y = probs)) + 
  geom_line(size=1) +
  xlim(-1.1, 0) +
  geom_vline(data = vogel_willems_results, 
             aes(xintercept = ate, linetype = study), 
             size=1,
             show.legend = F,
             color="grey50") + 
  geom_vline(xintercept = 0, size=1.5, color="firebrick") + 
  geom_text_repel(data=vogel_willems_results,
            aes(x = ate,
                y = seq(1, 0.75, length.out = 3),
                label = study),
                nudge_x = .05,
            segment.color = NA) + 
  labs(x = "Treatment effect on turnover intention",
       y = "Probability",
       title = "Prosocial recall task",
       linetype = "Study") + 
  scale_y_continuous(labels = percent_format()) + 
  theme_minimal_grid()) 

# Pulling out quantities of interest
less_than(variable = sim_eff_turnover$treatment_dikoProsocial,
             0)

# Probability of getting ATE of 0.3?
less_than(variable = sim_eff_turnover$treatment_dikoProsocial,
             -0.3)

# Probability of Vogel and Willems' ATEs?
tibble(
  study = vogel_willems_results$study,
  ate = vogel_willems_results$ate,
  probability =  str_c(100 * ate %>% 
  map_dbl(.f = less_than,
          variable = sim_eff_turnover$treatment_dikoProsocial), " pct."))

```

## Societal impact
```{r}
probs_societal <- seq_treat_size %>% 
  map_dbl(.f = less_than,
          variable = sim_eff_turnover$treatment_dikoSocietal)

results_societal <- cbind(seq_treat_size, probs_societal) %>% 
  as_tibble()

vogel_willems_results_societal <- 
  tribble(~study, ~ate,
          "Study 1", -1.05,
          "Study 2", -0.23,
          "Study 3", -0.68)

(p_prob_societal <- results_societal %>% 
  ggplot(aes(x = seq_treat_size, y = probs_societal)) + 
  geom_line(size=1) +
  xlim(c(-1.1, 0)) + 
  geom_vline(data = vogel_willems_results_societal, 
             aes(xintercept = ate, linetype = study), 
             size=1, 
             color="grey50",
             show.legend = F) + 
    geom_vline(xintercept = 0, size=1.5, color="firebrick") + 
  geom_text_repel(data=vogel_willems_results_societal,
            aes(x = ate,
                y = seq(1, 0.75, length.out = 3),
                label = study),
            nudge_x =  0.05,
            segment.color = NA) + 
  labs(x = "Treatment effect on turnover intention",
       y = "Probability",
       title = "Societal recall task") + 
  scale_y_continuous(labels = percent_format()) + 
  theme_minimal_grid()) 

# Probbability of any negative effect? (i.e. drop in turnover intention)
str_c(less_than(variable = sim_eff_turnover$treatment_dikoSocietal,
             0) * 100, " pct.")

# Probability of Vogel and Willems' ATEs?
tibble(
  study = vogel_willems_results_societal$study,
  ate = vogel_willems_results_societal$ate,
  probability =  str_c(100 * ate %>% 
  map_dbl(.f = less_than,
          variable = sim_eff_turnover$treatment_dikoSocietal), " pct."))
```

Putting the graphs together 
```{r}
full_plot_prob <- p_prob_prosoc / p_prob_societal

full_plot_prob + plot_annotation(tag_levels = 'A')

ggsave(filename= "prob_plot.tiff", 
       path="Figures",
       width = 10, 
       height = 10, 
       device='tiff', 
       dpi=700)
```

